/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

\*---------------------------------------------------------------------------*/

#include "immersedBoundaryFvPatchField.H"
#include "fvPatchFieldMapper.H"
#include "fvMatrix.H"
#include "surfaceWriter.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

template<class Type>
label immersedBoundaryFvPatchField<Type>::nBcIter_
(
    debug::optimisationSwitch("floatTransfernIter", 200)
);


template<class Type>
scalar immersedBoundaryFvPatchField<Type>::bcTolerance_
(
    1e-6//debug::tolerances("ImmersedBoundaryBCwTolerance", 1e-6)
);

// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

template<class Type>
void immersedBoundaryFvPatchField<Type>::updateIbValues() const
{
    // Create IB cell values first
    Field<Type> ibCellValues
    (
        this->internalField(),
        ibPatch_.ibCells()
    );

    if (this->fixesValue())
    {
        ibValue_ = ibPatch_.toIbPoints(refValue_);
        ibGrad_ = (ibValue_ - ibCellValues)/ibPatch_.ibDelta();
    }
    else
    {
        ibGrad_ = ibPatch_.toIbPoints(refGrad_);

        ibValue_ = ibCellValues + ibGrad_*ibPatch_.ibDelta();
    }
}


template<class Type>
Foam::tmp<Foam::Field<Type> >
immersedBoundaryFvPatchField<Type>::imposeDirichletCondition() const
{
    // Get addressing
    const labelList& ibc = ibPatch_.ibCells();
    const labelListList& ibCellCells = ibPatch_.ibCellCells();
    const List<List<labelPair> >& ibCellProcCells = ibPatch_.ibCellProcCells();
    const PtrList<scalarRectangularMatrix>& invMat =
        ibPatch_.invDirichletMatrices();

    const vectorField& ibp = ibPatch_.ibPoints();

    // Note: the algorithm is originally written with inward-facing normals
    // and subsequently changed: IB surface normals point outwards
    // HJ, 21/May/2012
    const vectorField& ibn = ibPatch_.ibNormals();

    // Collect Dirichlet values from IB triagulation
    ibValue_ = ibPatch_.toIbPoints(refValue_);
    
    // Reset the size and value of snGrad
    ibGrad_.setSize(ibc.size());
    ibGrad_ = pTraits<Type>::zero;

    // Get non-const access to internal field
    Field<Type>& psiI = const_cast<Field<Type>&>(this->internalField());

    // Collect polynomially interpolated values in IB cells
    tmp<Field<Type> > tpolyPsi(new Field<Type>(psiI, ibc));
    Field<Type>& polyPsi = tpolyPsi();

    const vectorField& C = mesh_.cellCentres();

    // Dimension the matrix
    label nCoeffs = 5;
    if (mesh_.nGeometricD() == 3)
    {
        nCoeffs += 4;
    }

    label counter = 0;
    scalarField error(ibc.size(), 0);

    do
    {
        counter++;

        // Parallel communication for psi
        FieldField<Field, Type> procPsi = ibPatch_.sendAndReceive(psiI);

        forAll (ibc, cellI)
        {
            label curCell = ibc[cellI];

            const labelList& curCells = ibCellCells[cellI];

            const List<labelPair>& curProcCells = ibCellProcCells[cellI];

            const scalarRectangularMatrix& curInvMatrix = invMat[cellI];

            Field<Type> coeffs(nCoeffs, pTraits<Type>::zero);
            Field<Type> source
            (
                curCells.size() + curProcCells.size(),
                pTraits<Type>::zero
            );

            label pointID = 0;
            for (label i = 0; i < curCells.size(); i++)
            {
                source[pointID++] = psiI[curCells[i]] - ibValue_[cellI];
            }

            for (label i = 0; i < curProcCells.size(); i++)
            {
                source[pointID++] =
                    procPsi
                    [
                        curProcCells[i].first()
                    ]
                    [
                        curProcCells[i].second()
                    ]
                  - ibValue_[cellI];
            }

            for (label i = 0; i < nCoeffs; i++)
            {
                for (label j = 0; j < source.size(); j++)
                {
                    coeffs[i] += curInvMatrix[i][j]*source[j];
                }
            }

            Type oldPolyPsi = polyPsi[cellI];

            vector R =  C[curCell] - ibp[cellI];

            polyPsi[cellI] =
                ibValue_[cellI]
              + coeffs[0]*R.x()
              + coeffs[1]*R.y()
              + coeffs[2]*R.x()*R.y()
              + coeffs[3]*sqr(R.x())
              + coeffs[4]*sqr(R.y());
            if (mesh_.nGeometricD() == 3)
            {
                polyPsi[cellI] +=
                    coeffs[5]*R.z()
                  + coeffs[6]*R.x()*R.z()
                  + coeffs[7]*R.y()*R.z()
                  + coeffs[8]*sqr(R.z());
            }
            
            // Change of sign of ibn
            ibGrad_[cellI] = (
               -coeffs[0]*ibn[cellI].x()
              - coeffs[1]*ibn[cellI].y());
            
            if (mesh_.nGeometricD() == 3)
            {
                // Change of sign of ibn
                ibGrad_[cellI] +=
                    -coeffs[5]*ibn[cellI].z();
            }

            error[cellI] = mag(polyPsi[cellI] - oldPolyPsi)
               /(mag(oldPolyPsi) + SMALL);
        }

        // Insert polynomial values into the internal field
        forAll (polyPsi, cellI)
        {
            psiI[ibc[cellI]] = polyPsi[cellI];
        }
    }
    while (gMax(error) > bcTolerance_ && counter < nBcIter_);

    if (counter == nBcIter_ && debug)
    {
        InfoIn
        (
            "template<class Type>\n"
            "tmp<Foam::Field<Type> >\n"
            "immersedBoundaryFvPatchField<Type>::"
            "imposeDirichletCondition() const"
        )   << this->dimensionedInternalField().name()
            << " for patch " << this->patch().name()
            << ", error, max: " << gMax(error)
            << ", min: " << gMin(error)
            << ", avg: "  << gAverage(error) << endl;
    }

//     Info<< "Dirichlet condition on " << ibc.size() << " cells of field "
//         << this->dimensionedInternalField().name() << " = "
//         << polyPsi
//         << endl;

    return tpolyPsi;
}


template<class Type>
Foam::tmp<Foam::Field<Type> >
immersedBoundaryFvPatchField<Type>::imposeNeumannCondition() const
{
    // Get addressing
    const labelList& ibc = ibPatch_.ibCells();
    const labelListList& ibCellCells = ibPatch_.ibCellCells();
    const List<List<labelPair> >& ibCellProcCells = ibPatch_.ibCellProcCells();
    const PtrList<scalarRectangularMatrix>& invMat =
        ibPatch_.invNeumannMatrices();

    const vectorField& ibp = ibPatch_.ibPoints();

    // Collect Neumann values from IB triagulation
    ibGrad_ = ibPatch_.toIbPoints(refGrad_);

    // Reset the size and value of
    ibValue_.setSize(ibc.size());
    ibValue_ = pTraits<Type>::zero;

    // Get non-const access to internal field
    Field<Type>& psiI = const_cast<Field<Type>&>(this->internalField());

    // Collect polynomially interpolated values in IB cells
    tmp<Field<Type> > tpolyPsi(new Field<Type>(psiI, ibc));
    Field<Type>& polyPsi = tpolyPsi();

    const vectorField& C = mesh_.cellCentres();

    // Dimension the matrix
    label nCoeffs = 6;
    if (mesh_.nGeometricD() == 3)
    {
        nCoeffs += 4;
    }

    label counter = 0;
    scalarField error(ibc.size(), 0);

    do
    {
        counter++;

        // Parallel communication for psi
        FieldField<Field, Type> procPsi = ibPatch_.sendAndReceive(psiI);

        forAll (ibc, cellI)
        {
            label curCell = ibc[cellI];

            const labelList& curCells = ibCellCells[cellI];
            const List<labelPair>& curProcCells = ibCellProcCells[cellI];

            const scalarRectangularMatrix& curInvMatrix = invMat[cellI];

            Field<Type> coeffs(nCoeffs, pTraits<Type>::zero);
            Field<Type> source
            (
                curCells.size() + 1 + curProcCells.size(),
                pTraits<Type>::zero
            );

            label pointID = 0;
            for (label i = 0; i < curCells.size(); i++)
            {
                source[pointID++] = psiI[curCells[i]];
            }

            source[pointID++] = ibGrad_[cellI];

            for (label i = 0; i < curProcCells.size(); i++)
            {
                source[pointID++] =
                    procPsi
                    [
                        curProcCells[i].first()
                    ]
                    [
                        curProcCells[i].second()
                    ];
            }

            for (label i = 0; i < nCoeffs; i++)
            {
                for (label j = 0; j < source.size(); j++)
                {
                    coeffs[i] += curInvMatrix[i][j]*source[j];
                }
            }

            Type oldPsi = polyPsi[cellI];

            vector ibR =  C[curCell] - ibp[cellI];

            polyPsi[cellI] = 
                coeffs[0]
              + coeffs[1]*ibR.x()
              + coeffs[2]*ibR.y()
              + coeffs[3]*ibR.x()*ibR.y()
              + coeffs[4]*sqr(ibR.x())
              + coeffs[5]*sqr(ibR.y());

            if (mesh_.nGeometricD() == 3)
            {
                polyPsi[cellI] +=
                    coeffs[6]*ibR.z()
                  + coeffs[7]*ibR.x()*ibR.z()
                  + coeffs[8]*ibR.y()*ibR.z()
                  + coeffs[9]*sqr(ibR.z());
            }
            ibValue_[cellI] = coeffs[0];

            error[cellI] = mag(polyPsi[cellI] - oldPsi)
               /(mag(oldPsi) + SMALL);
        }

        // Insert polynomial values into the internal field
        forAll (polyPsi, cellI)
        {
            psiI[ibc[cellI]] = polyPsi[cellI];
        }
    }
    while (gMax(error) > bcTolerance_ && counter < nBcIter_);

    if (counter == nBcIter_ && debug)
    {
        InfoIn
        (
            "template<class Type>\n"
            "tmp<Foam::Field<Type> >\n"
            "immersedBoundaryFvPatchField<Type>::"
            "imposeNeumannCondition() const"
        )   << this->dimensionedInternalField().name()
            << " for patch " << this->patch().name()
            << ", error, max: " << gMax(error)
            << ", min: " << gMin(error)
            << ", avg: "  << gAverage(error) << endl;
    }

//     Info<< "Neumann condition for " << ibc.size() << " cells of field "
//         << this->dimensionedInternalField().name() << " = "
//         << polyPsi
//         << endl;

    return tpolyPsi;
}


template<class Type>
void immersedBoundaryFvPatchField<Type>::imposeDeadCondition()
{
    const labelList& dc = ibPatch_.deadCells();

//     Info<< "Dead condition for " << dc.size()  << " cells of field "
//         << this->dimensionedInternalField().name()
//         << " set to value " << deadCellValue_
//         << endl;

    // Get non-const access to internal field
    Field<Type>& psiI = const_cast<Field<Type>&>(this->internalField());

    forAll (dc, dcI)
    {
        psiI[dc[dcI]] = deadCellValue_;
    }
}


template<class Type>
void immersedBoundaryFvPatchField<Type>::correctDiag
(
    fvMatrix<Type>& eqn
) const
{
    scalarField& Diag = eqn.diag();

    const labelList& dce = ibPatch_.deadCellsExt();

    // Estimate diagonal in live cells
    scalar liveDiag = 1;

    if (dce.size() < Diag.size())
    {
        liveDiag = gSumMag(Diag)/(Diag.size() - dce.size());

        // Correct for sign
        liveDiag *= sign(gMax(Diag));
    }

    forAll (dce, cellI)
    {
        if (mag(Diag[dce[cellI]]) < SMALL)
        {
            Diag[dce[cellI]] = liveDiag;
        }
    }
}


template<class Type>
void immersedBoundaryFvPatchField<Type>::correctOffDiag
(
    fvMatrix<Type>& eqn
) const
{
    // Calculate gradient contribution
    const labelList& ibFaces = ibPatch_.ibFaces();
    const labelList& ibFaceCells = ibPatch_.ibFaceCells();

    const scalarField& ibGamma = ibPatch_.gamma().internalField();

    const unallocLabelList& own = mesh_.owner();
    const unallocLabelList& nei = mesh_.neighbour();

    // Get delta coefficients
    const surfaceScalarField& dc = mesh_.deltaCoeffs();
    const scalarField& dcI = dc.internalField();

    if (eqn.symmetric())
    {
        scalarField& diag = eqn.diag();
        scalarField& upper = eqn.upper();
        Field<Type>& source = eqn.source();

//         Info<< "Symmetric correctOffDiag for field "
//             << this->dimensionedInternalField().name() << endl;

        forAll (ibFaces, faceI)
        {
            const label curFace = ibFaces[faceI];

            if (curFace < nei.size())
            {
                // Internal face.  One side is an ibCell and another is a
                // live cell. Add gradient to the source of the live cell
                // and kill the off-diagonal coefficient
                if (ibGamma[own[curFace]] > SMALL)
                {
                    diag[own[curFace]] += upper[curFace];

                    source[own[curFace]] +=
                        upper[curFace]*ibGrad_[ibFaceCells[faceI]]
                        /dcI[curFace];
                }
                else
                {
                    diag[nei[curFace]] += upper[curFace];

                    source[nei[curFace]] -=
                        upper[curFace]*ibGrad_[ibFaceCells[faceI]]
                        /dcI[curFace];
                }

                upper[curFace] = 0;
            }
            // else MISSING PARALLELISATION HANDLING
        }
    }
    else if (eqn.asymmetric())
    {
        scalarField& diag = eqn.diag();
        scalarField& upper = eqn.upper();
        scalarField& lower = eqn.lower();
        Field<Type>& source = eqn.source();

//         Info<< "Asymmetric correctOffDiag for field "
//             << this->dimensionedInternalField().name() << endl;

        forAll (ibFaces, faceI)
        {
            const label curFace = ibFaces[faceI];

            if (curFace < nei.size())
            {
                // Internal face.  One side is an ibCell and another is a
                // live cell. Add gradient to the source of the live cell
                // and kill the off-diagonal coefficient
                if (ibGamma[own[curFace]] > SMALL)
                {
                    diag[own[curFace]] += upper[curFace];

                    source[own[curFace]] +=
                        upper[curFace]*ibGrad_[ibFaceCells[faceI]]/dcI[faceI];
                }
                else
                {
                    diag[nei[curFace]] += lower[curFace];

                    source[nei[curFace]] -=
                        lower[curFace]*ibGrad_[ibFaceCells[faceI]]/dcI[faceI];
                }

                upper[curFace] = 0;
                lower[curFace] = 0;
            }
            // else MISSING PARALLELISATION HANDLING
        }
    }


    // Note: potentially deal with face flux correction ptr.
    // HJ, 16/Apr/2012
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class Type>
immersedBoundaryFvPatchField<Type>::immersedBoundaryFvPatchField
(
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF
)
:
    fvPatchField<Type>(p, iF, Field<Type>(0)),
    ibPatch_(refCast<const immersedBoundaryFvPatch>(p)),
    mesh_(p.boundaryMesh().mesh()),
    refValue_(ibPatch_.ibMesh().size(), pTraits<Type>::zero),
    refGrad_(ibPatch_.ibMesh().size(), pTraits<Type>::zero),
    fixesValue_(false),
    setDeadCellValue_(false),
    deadCellValue_(pTraits<Type>::zero),
    ibValue_(),
    ibGrad_()
{}


template<class Type>
immersedBoundaryFvPatchField<Type>::immersedBoundaryFvPatchField
(
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF,
    const dictionary& dict
)
:
    fvPatchField<Type>(p, iF, Field<Type>(0)),
    ibPatch_(refCast<const immersedBoundaryFvPatch>(p)),
    mesh_(p.boundaryMesh().mesh()),
    refValue_("refValue", dict, ibPatch_.ibMesh().size()),
    refGrad_("refGradient", dict, ibPatch_.ibMesh().size()),
    fixesValue_(dict.lookup("fixesValue")),
    setDeadCellValue_(dict.lookup("setDeadCellValue")),
    deadCellValue_(pTraits<Type>(dict.lookup("deadCellValue"))),
    ibValue_(),
    ibGrad_()
{
    if (!isType<immersedBoundaryFvPatch>(p))
    {
        FatalIOErrorIn
        (
            "immersedBoundaryFvPatchField<Type>::"
            "immersedBoundaryFvPatchField\n"
            "(\n"
            "    const fvPatch& p,\n"
            "    const Field<Type>& field,\n"
            "    const dictionary& dict\n"
            ")\n",
            dict
        )   << "\n    patch type '" << p.type()
            << "' not constraint type '" << typeName << "'"
            << "\n    for patch " << p.name()
            << " of field " << this->dimensionedInternalField().name()
            << " in file " << this->dimensionedInternalField().objectPath()
            << exit(FatalIOError);
    }
}


template<class Type>
immersedBoundaryFvPatchField<Type>::immersedBoundaryFvPatchField
(
    const immersedBoundaryFvPatchField<Type>& ptf,
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF,
    const fvPatchFieldMapper&
)
:
    fvPatchField<Type>(p, iF, Field<Type>(0)),
    ibPatch_(refCast<const immersedBoundaryFvPatch>(p)),
    mesh_(p.boundaryMesh().mesh()),
    refValue_(ptf.refValue()),
    refGrad_(ptf.refGrad()),
    fixesValue_(ptf.fixesValue()),
    setDeadCellValue_(ptf.setDeadCellValue_),
    deadCellValue_(ptf.deadCellValue_),
    ibValue_(),
    ibGrad_()
{
    // Note: NO MAPPING.  Fields are created on the immersed boundary
    // HJ, 12/Apr/2012
    if (!isType<immersedBoundaryFvPatch>(p))
    {
        FatalErrorIn
        (
            "immersedBoundaryFvPatchField<Type>::"
            "immersedBoundaryFvPatchField\n"
            "(\n"
            "    const immersedBoundaryFvPatchField<Type>&,\n"
            "    const fvPatch& p,\n"
            "    const DimensionedField<Type, volMesh>& iF,\n"
            "    const fvPatchFieldMapper& mapper\n"
            ")\n"
        )   << "\n    patch type '" << p.type()
            << "' not constraint type '" << typeName << "'"
            << "\n    for patch " << p.name()
            << " of field " << this->dimensionedInternalField().name()
            << " in file " << this->dimensionedInternalField().objectPath()
            << exit(FatalIOError);
    }
}


template<class Type>
immersedBoundaryFvPatchField<Type>::immersedBoundaryFvPatchField
(
    const immersedBoundaryFvPatchField<Type>& ptf
)
:
    fvPatchField<Type>
    (
        ptf.patch(),
        ptf.dimensionedInternalField(),
        Field<Type>(0)
    ),
    ibPatch_(ptf.ibPatch()),
    mesh_(ptf.patch().boundaryMesh().mesh()),
    refValue_(ptf.refValue()),
    refGrad_(ptf.refGrad()),
    fixesValue_(ptf.fixesValue()),
    setDeadCellValue_(ptf.setDeadCellValue_),
    deadCellValue_(ptf.deadCellValue_),
    ibValue_(),
    ibGrad_()
{}


template<class Type>
immersedBoundaryFvPatchField<Type>::immersedBoundaryFvPatchField
(
    const immersedBoundaryFvPatchField<Type>& ptf,
    const DimensionedField<Type, volMesh>& iF
)
:
    fvPatchField<Type>(ptf.patch(), iF, Field<Type>(0)),
    ibPatch_(ptf.ibPatch()),
    mesh_(ptf.patch().boundaryMesh().mesh()),
    refValue_(ptf.refValue()),
    refGrad_(ptf.refGrad()),
    fixesValue_(ptf.fixesValue()),
    setDeadCellValue_(ptf.setDeadCellValue_),
    deadCellValue_(ptf.deadCellValue_),
    ibValue_(),
    ibGrad_()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class Type>
const Field<Type>& immersedBoundaryFvPatchField<Type>::ibValue() const
{
    if (ibValue_.empty())
    {
        this->updateIbValues();
    }

    return ibValue_;
}


template<class Type>
const Field<Type>& immersedBoundaryFvPatchField<Type>::ibGrad() const
{
    if (ibGrad_.empty())
    {
        this->updateIbValues();
    }

    return ibGrad_;
}

template<class Type>
tmp<Field<Type> > immersedBoundaryFvPatchField<Type>::triValue() const
{
    return ibPatch_.toTriFaces(this->ibValue());
}


template<class Type>
tmp<Field<Type> > immersedBoundaryFvPatchField<Type>::triGrad() const
{
    return ibPatch_.toTriFaces(this->ibGrad());
}


template<class Type>
void immersedBoundaryFvPatchField<Type>::updateCoeffs()
{
    if (this->fixesValue())
    {
        this->imposeDirichletCondition();
    }
    else
    {
        this->imposeNeumannCondition();
    }

    // Fix the value in dead cells
    if (setDeadCellValue_)
    {
        this->imposeDeadCondition();
    }

    fvPatchField<Type>::updateCoeffs();
}


template<class Type>
void immersedBoundaryFvPatchField<Type>::initEvaluate
(
    const Pstream::commsTypes
)
{}


template<class Type>
void immersedBoundaryFvPatchField<Type>::evaluate
(
    const Pstream::commsTypes
)
{
    // Note
    // Since the boundary condition is performed by data fitting with the
    // internal field, fitting must be performed both on updateCoeffs
    // and on evaluate (internal field has changed in the meantime).
    // Bug fix, Zeljko Tukovic, 21/Jun/2012
    this->updateCoeffs();

    fvPatchField<Type>::evaluate();
}


template<class Type>
void immersedBoundaryFvPatchField<Type>::manipulateMatrix
(
    fvMatrix<Type>& eqn
)
{
    this->initEvaluate();

    // Build matrix diagonal for cells where it is missing
    this->correctDiag(eqn);

    // For Neumann boundary condition, manipulate matrix off-diagonal
    if (!this->fixesValue())
    {
        this->correctOffDiag(eqn);
    }

    // Set values in IB cells
    Field<Type> polyPsi(eqn.psi(), ibPatch_.ibCells());
    eqn.setValues(ibPatch_.ibCells(), polyPsi);

    // Correct equation for dead cells
    Field<Type> deadCellsPsi
    (
        ibPatch_.deadCells().size(),
        deadCellValue_
    );
    eqn.setValues(ibPatch_.deadCells(), deadCellsPsi);

    fvPatchField<Type>::manipulateMatrix(eqn);
}


template<class Type>
void immersedBoundaryFvPatchField<Type>::write(Ostream& os) const
{
    fvPatchField<Type>::write(os);
    //agregado
    cout<<"refValue_ MODIFICADO para enganchar con dynamicBodyRefineFvMesh - NO VALIDO PARA OTROS CAMPOS!\n";
    Field<Type>& rvAux = const_cast<Field<Type>&>(refValue());
    Field<Type> rvOrig(rvAux);
    rvAux *= 0;
    refValue_.writeEntry("refValue", os);
    rvAux = rvOrig;
    
    refGrad_.writeEntry("refGradient", os);
    os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
    os.writeKeyword("setDeadCellValue")
        << setDeadCellValue_ << token::END_STATEMENT << nl;
    os.writeKeyword("deadCellValue")
        << deadCellValue_ << token::END_STATEMENT << nl;

    this->writeEntry("value", os);

    // Write immersed boundary data as a vtk file
    /*autoPtr<surfaceWriter<Type> > writerPtr =
        surfaceWriter<Type>::New("vtk");

    const triSurface& ts = ibPatch_.ibMesh();

    // Make a face list for writing
    faceList f(ts.size());
    forAll (ts, faceI)
    {
        f[faceI] = ts[faceI].triFaceFace();
    }

    writerPtr->write
    (
        this->dimensionedInternalField().path(),
        ibPatch_.name(),
        ts.points(),
        f,
        this->dimensionedInternalField().name(),
        this->triValue()()
    );*/
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //
